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Abstract: Survey sampling methods provide cost-effective solutions for monitoring 
levels in large populations. Rotation designs, which replace a fraction of the sam- 
ple at each survey occasion, yield better performances than fixed panels for this 
purpose. Motivated by recent applications in sensor networks, we extend rotation 
sampling to functional data (that is, continuous signals). We propose survey de- 
signs that replace the sample at fixed times via Markov chains. We then study the 
related Horvitz- Thompson estimator of the population mean function and quantify 
its large-sample behavior. For a given sample size, the mean estimation error can be 
decreased by suitably allocating the sample during replacements; the error variance 
can be dramatically reduced by increasing either the frequency or the intensity of 
replacements. We also extend the composite estimation method to the continuous- 
time setup. The new estimator combines the Horvitz-Thompson estimator, which 
only uses current data, with the estimated population change between an arbitrary 
past time and the current observation time. A simulation study shows that the pro- 
posed rotation designs and composite estimator considerably improve upon fixed 
panels and the Horvitz-Thompson estimator. 

Key words and phrases: Design-based survey; functional data; rotation sampling; 
Horvitz-Thompson estimator; asymptotic theory; composite estimator. 

1. Introduction 

In various industrial, environmental and medical applications, sensor net- 
works generate large volumes of data in continuous time. Due to cost or energy 
constraints, such collections of functional data (that is, curve data) often cannot 
be entirely observed. For instance, large electric utilities need to monitor the 
total consumption of their clients in order to adjust the power generation to the 
system load, predict future consumption, and determine pricing policies. These 
utilities however cannot read all of their client smart meters at each instant, as 
the network cannot process such large data transmission and the cost for data 
storage would be prohibitive. Under such observation constraints, survey sam- 
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pling methods offer competitive solutions for monitoring global parameters; see 
e.g., Chiky, Cubille, Dessertaine, Hebrail, and Picard (2008) for a comparison 
between survey sampling and signal compression approaches. 

Repeated surveys are well studied in the statistical literature. In the clas- 
sical modeling framework where population parameters are fixed and only the 
sample selection is random, authors like Patterson (1950), Eckler (1955), Rao 
and Graham (1964), and Wolter (1979) have developed the theory of rotation 
sampling, a popular technique that consists in replacing a fraction of the sample 
with new population units at each survey occasion. They have established that 
different replacement strategies should be used depending on the survey's goal. 
For example, partial replacements of the sample are advantageous for assessing 
population levels while fixed panels are preferable for estimating level changes. 
Jointly with the study of sampling designs, they investigated estimation proce- 
dures such as minimum variance unbiased estimators and composite estimators, 
aiming to strike a balance between statistical and computational efficiency. 

While traditional repeated surveys take place in discrete time and present 
important limitations on the sampling design and data collection, sensor networks 
such as smart electrical grids or audience measurement systems provide data ev- 
ery few seconds or minutes (essentially, in continuous time). This recent type of 
application enables sophisticated sampling designs and alleviates burdens of tra- 
ditional surveys such as measurement errors, nonresponse and sample attrition. 
At the same time, it poses new statistical challenges related to the functional 
data/continuous-time setup. For example, Cardot and Josserand (2011) and 
Cardot, Degras, and Josserand (2012) extend the well-known Horvitz-Thompson 
[HT] estimator to continuous-time surveys. Based on results from Degras (2011) 
they address the construction of simultaneous confidence bands, which of course 
only makes sense in view of a continuous data-generating process. Studying the 
model-assisted survey of functional data, Cardot, Chaouch, Goga, and Labruere 
(2010) use functional principal component analysis and multivariate nonparamet- 
ric regression to relate scalar auxiliary information to functional observations. 
Note that the above survey methods are all based on fixed panels. 

In this paper we devise rotation designs for the survey of functional data. 
In contrast with traditional rotation designs which specify the sample before the 
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survey (see e.g., Eckler, 1955), we build rotation samples using Markov chains. 
As in Tikkiwal and Gupta (1991) this approach introduces a new randomization 
at each sample replacement. Importantly, it enables informative sampling: past 
measurements can be used to "improve" the sample at the next replacement by 
suitably updating the inclusion/exclusion probabilities. Examples of adaptive 
sampling rules would be to retain units in the sample with probabilities propor- 
tional to their current levels, or to optimize the sample allocation in stratified 
sampling based on the current strata sizes and variances. 

In a second time, we examine the estimation of the population mean function 
based on our proposed rotation designs. Specifically, we quantify the large-sample 
behavior of the HT estimator in terms of the sampling rate, the replacement rate 
(that is, the fraction of the sample discarded at each replacement), and the 
number of replacements. We focus on the integrated squared error [ISE] which 
is a relevant measure for assessing a sampling strategy over a time interval. Our 
results show that when the sample size is constant over time, rotation samples 
and fixed panels yield the same mean ISE. On the other hand, the mean ISE can 
be sensibly reduced under rotation sampling by using stratification and suitably 
re-allocating the sample at replacement times. More importantly, the variance 
of the ISE can be decreased by several orders of magnitude when using rotation 
samples rather than fixed panels. This decrease is achieved by increasing either 
the frequency or the intensity of replacements. 

Thirdly, this paper extends the composite estimation method (see e.g., Rao 
and Graham, 1964) to functional data. Unlike discrete-time surveys where com- 
posite estimators are updated with respect to the last survey occasion, there is 
no "last" occasion in continuous time. We thus define our composite estimator 
in reference to an arbitrary past observation time. This estimator is based on 
the HT estimator and an estimator of the population change. Note that the the- 
ory of minimum variance unbiased estimation (see e.g., Patterson, 1950) is of no 
help in our context as we do not consider infinite populations nor make assump- 
tions on the temporal dependence. To finish, we present a numerical study with 
artificial electricity consumption data that confirms the superior performances 
of our rotation samples and composite estimator over fixed panels and the HT 
estimator. 
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The rest of the paper is organized as follows. The modeling framework 
and the HT estimator are presented in Section 2. The new rotation designs for 
functional data are defined in Section 3. The mean and covariance of the HT 
estimator are studied in Section 4. Section 5 contains our main result, namely 
the large-sample variance of the ISE under the studied rotation designs. The 
composite estimator is introduced in Section 6 and the simulation study is de- 
scribed in Section 7. Concluding remarks are offered in Section 8. The proofs 
of the main results are collected in the Appendix and other proofs are available 
online as supplementary material. 

2. Statistical framework 

Let Un = {1, • • • ,N} be a finite population in which a deterministic curve 
Xk(t),t G [0, T], is associated to each unit k G Un- We study the estimation of 
the population mean function 

keU N 

based on observations -Xfc(i), k G s(t), where s(t) C Un is a probability sample of 
fixed size n(t). For mathematical simplicity, we assume perfect observations. If 
the sampled curves are observed at discrete times and/or with noise, interpolation 
or smoothing methods should be applied to the data. In this case the results of 
this paper will still hold under standard interpolation or nonparametric smooth- 
ing conditions. See for example Cardot and Josserand (2011) and Cardot, Degras 
and Josserand (2012). The sample s(-) = {s(t), t G [0,T]} is selected according 
to a controlled probability distribution pn over the function space V(Un)^°' T ^ , 
where V(Un) denotes the set of all subsets of Un- For any units k,l G Un and 
times t,t' G [0, T], the first and second order inclusion probabilities under pn are 
written as ir^it) = G s(t)) and TTki(t,t') = F(k £ s(t), I G s(t')) respectively. 
Henceforth the subscript in Un is dropped. 

We use the estimator of Horvitz and Thompson (1952) 

where Ik(t) is the sample indicator function of k G Un at time t: Ik(t) = 1 if 
k G s(t) and /&(£) = otherwise. The HT estimator is unbiased under pn and 
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its covariance function is 

Cov(A«(t),A«(0) = ^$) Xk{t)Xl{t% 

where A w (*,f) = Cov(/ fc (t),/,(f )) = 7r fc ,(t,0 - 7r fc (t)ir,(/ / ). 

The estimation accuracy can be improved by using a suitable stratification 
of the population. From now on we assume that U is partitioned into strata 
of size Nh for h = 1, . . . , if. For each £7^, denote by //^(t) = (1/Nh) ^2keu h -^k(t) 
and lh {t, t') = {l/(N h - 1)) J2keu h ( X k(t) ~ A*h(*)) (*k(0 " Mfc(O) th e stratum 
mean and covariance functions. Let n^it) = #(s(t) D Uh) be the sample size in 
Uh at time t and fh(t) = rih(t)/Nh be the sampling rate. If s(t) is obtained by 
simple random sampling without replacement [SRSWOR] independently in each 
Uh, the HT estimator becomes 

1 H 1 

Mt) = - J7T) E h(t)X k (t) (2.2) 

h=l Jh{ - ' keU h 

and its covariance rewrites as 

1 H 1 

cov(A ftt (t),Ahf(0) = f (t )f m E A„(i, 0^(^(0- (2-3) 

3. Rotation designs for continuous-time surveys 

The rotation sampling technique has been developed in the context of (discrete- 
time) traditional repeated surveys. In this section we extend rotation sampling to 
the continuous-time setup of functional data. We propose two novel probability 
schemes for selecting a time- varying sample s(-) = {s(t),t G [0, T]} in a stratified 
population. These sampling designs, which we refer to as full replacement and 
partial replacement, share the following features: 

• The strata samples Sh(-) = {s(t) C\ Uh,t G [0, T]} are independent across 
strata. 

• At time to = 0, the strata samples Sfc(ro) are obtained by SRSWOR. 

• The Sh( ) can only be modified at fixed times < t\ < . . . < r m < T. 
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It remains to specify the evolution mechanisms of the discrete processes {s/ l (r r ), 
r = 1, . . . , m} under full and partial replacement. 

1. Full replacement. In each stratum Uh, the successive samples s^(r r ), r = 
1, . . . , m, are obtained by independent SRSWOR of n/ l (r r ) units in Uh- 

2. Partial replacement. In each stratum Uh, a fraction ah £ [0,1] of Sh(-) 
is replaced at each time r r , r = 1, . . . , m. Specifically, s^(r r _i) is updated 
at time r r by the following independent operations: 

• discard ahrih(T r -i) units selected in s^(r r _i) by SRSWOR; 

• add (n h (T r ) - n/ l (r r _i) + a h n h (T r -i)) units selected in U h \ s h {T r -i) by 
SRSWOR. 

Note that full replacement is not a special case of partial replacement with ah = 1. 
Indeed, for r = 1, . . . , m, Sh{T r -%) and Sh(r r ) are independent in the former case 
whereas they are disjoint (and thus dependent) in the latter case. In partial 
replacement we refer to the ah as the replacement rates. For simplicity we as- 
sume that the ah are constant over time and that the proposed sample replace- 
ments are possible without modifications, which entails that ahnh( T r-i) £ N and 
n^(r r _i) < nh(r r ) + Q/ i n/ 1 (r r _i) < Nh for all h,r. Note that fixed panels are a 
special case of partial replacement with ah = and rih(-) constant over time. 

We now determine the probability distribution of the sample s(t) under the 
proposed sampling designs. The following result relies on an induction argument 
on the r r under partial replacement while it holds trivially under full replacement. 

Proposition 1. Consider either the full or the partial replacement design. For 
each stratum Uh and time t S [0, T], the probability distribution of the sample 
Sh(t) is identical to the SRSWOR of nh{t) units in Uh- 

4. Covariance of the Horvitz-Thompson estimator 



Here we derive the covariance function of the HT estimator (2.2) under the 
full and partial replacement designs of the previous section. Let v{t) denote the 
number of sample replacements occurring before time t: 

u(t) = min {r G {0, . . . , m} : r r < t} 

so that t £ [tl,(£\, for all t G [0, T) with the convention r m+ \ = T. 
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4.1. Covariance under the full replacement design 

Under the full replacement design, fihtif) an d Aht(*0 are independent if the 
sample has been replaced between times t and t' . Therefore the expression of the 



covariance (2.3) directly follows from the properties of SRSWOR. Let 5.. indicate 



the Kronecker delta. 

Theorem 1. Consider the full replacement design of Section 3. For all strata 
Ufi, units k,l E Uh, and times t,t' S [0, T], 

A«(t, = (1 - f h {t)) f h (t) K {t)u{t>) iY ^~ 1 • 



As a consequence, the covariance of the HT estimator (2.2) is 

Cov(/*w(*),Am(0) = ]y ^ -Jf JJ^j 7h(*. t') 5 v ( t ) V (t>) ■ 

This theorem will be commented in relation to partial replacement in the next 
section. 

4.2. Covariance under the partial replacement design 

In view of Proposition [TJ it suffices to compute 7r^(t, t') = E(ifc(t)/;(i')) for 
k,l £ Uh (h = 1, . . . , H) and < t < t' < T in order to determine Afc/(i, t') and 



(2.3). By definition of SRSWOR, the inclusion probabilities 7Tfcfc(i, i') and 7Tfcz(t, t') 
do not depend on G f7/j (fc 7^ /). Since 5Zfc-^fc(*) = n hit), one also sees that 
Cov (E fce c4 Ikit),EieU h W)) = N h A kk it,t') + N h iN h -l) Ay(M') = 0, where 
7^ / are two arbitrary units in Uh- Therefore, computing Ajy(t, t') amounts to 
computing A^(t,t') which in turn reduces to deriving P(fe G s/j(i')|fc E Shit))- 

The transition probabilities of the Markov chain {/fc(r r ), r = 0, . . . ,m} can 
be found by applying the Chapman-Kolmogorov equations. To that end, define 

MM')= n ; ah ~ f w 

1 - fh{Tr-l) 

for all times < t < t' < T and extend A/ l (t, i') as a symmetric function on 
[0,T] 2 . Set A/ l (t,f) = 1 if u(t) = u(i?) and set to 1 all factors of X h it,t') for 
which /h(r r _i) = 1. 
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Lemma 1. Consider the partial replacement design of Section 3. For all strata 
Uh, units k £ Uh, and times < t < t' < T , 

f P (k E s h (t')\k e s h (t)) =(l-f h (t))X h (t,t') + f h (t'), 

\ p (k e s h (t')\k i Sh (t)) = f h (t') - f h (t)\ h (t,t>). 

In the proof of this lemma, the quantity Xh(t, t') turns out to be the product 
of the eigenvalues of the transition probability matrices of the previous Markov 
chain between times t and t'. 

For any two real numbers x, y, write xAy = min(x, y) and x\/y = max(x, y). 
Exploiting Proposition [T] and Lemma [TJ we obtain the covariance function (2.3) 
under partial replacement. 

Theorem 2. Consider the partial replacement design of Section 3. For all strata 
Uh, units k,l E Uh, and times t,t' £ [0, T], 

A w (t, = (1 - f h (t A t')) f h {t A jV ^~ 1 X h (t, t'). 



As a consequence, the covariance of the HT estimator (2.2) is 

Cov(A«(t),A«(0) = Jf E § 1 ; ft {"v t'P lh{t ' * Xhit ' t0 " 

To provide insights into Theorems T][2 we consider the situation where the 
sample sizes n/ l (-) are constant over time. In this case the term \h(t,t') in 
Theorem [i] simplifies into (1 — a/t/(l — A))' 1 '^ -1 ^* We can then compare the 
estimator covariance under the full and partial replacement designs, examining 
the "no replacement" strategy (that is, a fixed panel) as a special instance of 
partial replacement (ah = 0). On the diagonal blocks [r r , r r+ i] 2 , r = l,...,m, 
Cov(£iht(t), Aht(i')) is identical under the three previous types of replacement. 
Outside these blocks, the estimator covariance is zero under full replacement 
whereas it is a weighted average of the covariance functions 7/1 under a fixed panel. 
For partial replacement the covariance structure depends on the replacement 
rates ah- When ah increases in [0, 1 — fh] for some h, the estimator covariance 
between given times t, t' decreases (assuming positivity of 7/J. In the special case 
where ah = 1 — fh for all h, the covariance is the same as under full replacement. 
For values ah > 1 — fh, the covariance becomes unstable in the sense that the 
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term (1 - a h j (1 - jh]j v{ ^~ v ^ ^ changes sign on every block [r q , T q+1 ] x [r r , r r+ i]. 
For ah $l {0, 1 — fh} and regularly spaced r r (that is, \v{t) — v{t')\ < Cm\t — t'\ 
for some finite constant C), the estimator covariance decreases at an exponential 
rate as \t — t'\ increases. 

4.3. Mean Integrated Squared Error 

To measure the accuracy of an estimator ji^ of fijy over [0, T], we use the 
integrated squared error 

ISE= / (fi N (t) - m (t)) 2 dt. 



As seen in Sections 2-3, the HT estimator (2.2) is unbiased and, when the sample 
sizes rih(-) are constant over time, its variance function is the same under the full 
and partial replacement designs (and in particular under fixed panels). In this 
case, the HT estimator has the same mean integrated squared error 

rT 



MISE= f E(fi N (t)- ^i N {t)) 2 dt 
Jo 



under both designs. The MISE can however be reduced (in comparison to fixed 
panels) by using the full or partial replacement designs with suitable time- varying 
sample sizes n^(-). Specifically, the variance of fiht(t) c an be minimized at 
each replacement time r r by choosing n/ l (T r ) according to the classical Ney- 
man allocation rule n/i(r r ) = c r Nh \f"1h{ T ri T r) with the constant c r such that 
^2h n h( T r) = n(r r ). See for example Fuller [2009, p. 21] for more details. 

5. Asymptotic results for the ISE 



We now determine the variance of the ISE for the HT estimator (2.2) based 
on the full or partial replacement designs of Section 3. Starting with the general 



HT estimator (2.1), this variance expresses as 



Var (ISE) = JJ^ Cov ({£**(*) - ^ N (t)} 2 , {fi ht {t') - ^ N (t')} 2 ) dtdt' 

1 ff V ^IMMM) A (tt')dtdt' (^) 



where we have introduced the four-fold cross-covariance function 



&ijkl{t, t') = Cov({Ii(t) - ni(t)}{Ij(t) - TTj (t)} , {h(t') - KktfMW) - Tl(0})- 
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While (5.1) can be computed exactly if the sample sizes rih(-) are constant 
over time, large-sample approximations are required under time- varying n/ l (-). 
The asymptotic framework of such approximations is detailed in Section 5.1, 
intermediate results are presented in Section 5.2, and the main results are given 
in Section 5.3. 

5.1. Asymptotic framework 

We let the strata sizes Nh, sample sizes n^(-), replacement rates ah (h = 
1, . . . , H), and number of replacements m depend on the population size iV and 
go to infinity (except for ah) together with N. The number H of strata and the 
observation period [0, T] stay fixed. We make the following assumptions. 

(Al) The curves Xk, k>l, are integrable and uniformly bounded on [0, T]. 

I-Tr 

(A2) / g(t)dt = r/(m + 1) for r = 1, . . . , m, where g is a continuous density 
Jo 

function supported by [0, T]. 

(A3) The sampling rate functions fh(-) stay bounded away from zero as N — > oo. 

(A4) The covariance functions 7^ converge uniformly on [0,T] 2 as N — > 00 and 
have continuous limits. 

(A5) The number of replacements is dominated by the strata sizes: m = o(min/ l (A^/ 1 )) 
as N — > 00. 

The number H of strata, although fixed, can be large. The condition Nh — > 00 
is not restrictive as, typically, small strata Uh are fully observed and do not 
contribute to the estimation error. In (Al) the individual curves are allowed 
to have discontinuity jumps. However (A4) requires that the strata covariance 
functions can be uniformly approximated by continuous functions, which entails 
that at any time t, at most an asymptotically negligible fraction of the X^{t) may 
have discontinuity jumps. This assumption is needed under full replacement 
to approximate 7h(t, t') around the diagonal {t = t'} by the variance jh{t,t). 
Assumption (A2) ensures that the replacement times are regularly spaced, while 
(A3) is necessary for the consistent estimation of //jy. Finally (A5) is needed 
under partial replacement to approximate the transition probabilities in and out 
of the sample for fixed population units. 
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5.2. Intermediate results 

Here we derive the cross-covariance function Aijki(t,t'). Recall that the 
time- varying samples Sh(-) are independent across strata. By expressing the 
population mean function fijy as the weighted average Ylhi^h/^l^h of the strata 
mean functions and using the corresponding decomposition for the HT estimator 
( [272] ), it follows that 

(5.2) 

Writing Xk(t) = Xk(t) — /J,h(t) for k G Uh and h = 1, . . . , H, simple algebra 
shows that 

A -« (t, X t (t)Xj (t)X k (t')Xi (t') 

i,j,k,leU h 

= Yl E(/ i (t)J i (t)/ fc (t')J,(0)X i (t)X i (*)X Jfe (t / )^i(0 ( 5 - 3 ) 

i,j,k,leUh 

- N 2 h f h (t)f h (t>) (1 - f h {t)) (1 - A(t')) 7h(*, t) J h (t', t') . 



The sum in the right-hand side of (5.3) can be developed using the properties 
of SRSWOR. We denote the asymptotic equivalence of two real sequences (ajv) 
and (&at) by ~ b^. 

Proposition 2. Consider either the full or the partial replacement design of 
Section 3. Assume (Al). Let {i* , j* , k* ,1*} be four distinct units in a given 
stratum Uh- Then 

i,j,k,laUh 

~ (Ci(t, 7 fc(t, 07ft (*', + C 2 (t, t') 7 2 h (t, t')) N 2 h 
uniformly in t, t' G [0, T] as iV^ —> oo, where 

d(t,0 = E (lHt)4*(0) - E (li,(t)I r (t)I k ,(t')) - E (/ i .(t)/ft.(t')I l .(* / )) 
+ E(J i .(*)J i .(*)J fc .(t')J,.(0) 

and 

C 2 (M') = 2 E (If. (t)ii. (i)Jfc. (i')) -4E (Ji.(t)Ji. (t')Ij*(t)I k * (tO) 
+ 2E (/<.(*)/,-.(«) Jfc.Ct')/,.^)). 
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Under the full replacement design, the previous functions G\ and C2 can be 
expressed in terms of fh(t) and fh(t') and Var(ISE) can readily be computed. 
For partial replacement however, additional results are required. 

We now provide a result on the distribution of the sample in a given subset 
conditional on past values of the sample. Denote conditional distributions by 
£(-|-). Fix a stratum Uh and a subset D C Uh- The Markovian nature of 
{sh( T r), f = 0, . . . , m} and the properties of SRSWOR (namely, the probability 
that the sample contains D only depends on the size of D) yield the following 
result. 



Lemma 2. Consider the partial replacement design of Section 3. Then for all 
<t <t' <T, 



The lemma states that at time t, all information relative to the future distribution 
of Sh(t') n D is contained in Sh(t) n D. 

Using the Chapman-Kolmogorov equations, large-sample approximations and 
linear algebra, we derive the transition probabilities in and out of the sample for 
two units of a given stratum. 

Proposition 3. Consider the partial replacement design of Section 3 and assume 
(A3)-(A5). For all units k,l £ Uh (k ^ I) and times <t < t' <T, it holds as 
N — > 00 that 



5.4. Variance of the Integrated Squared Error 

Based on the previous findings, we can now state the main results. 

Theorem 3. Consider the HT estimator based on the full replacement 

design of Section 3. Assume (A1)-(A2)-(A3)-(A4). Then as N -4 00, 



c(s h (t') n d I Sh (t)) = c(s h (t') n d I s h (t) n d). 



' P(k,l€ s h (t')\k,le s h (t)) ~ [(l-f h (t))X h (t,t') + f h (t')] 2 , 
P (k,l 6 s h (t')\k e s h (t), I $ s h {t)) 



< 



~ [ - f h (t) (1 - f h {t)) \l(t, t') + f h (t>) (1 - 2 f h (t)) \ h it) + f h it')] 

k P(M € «h(0|M £ */»(«)) ~ [(1- hit)) \ h (t,t') - (1 - h(t'))] 2 . 
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Theorem 4. Consider the HT estimator (2.2) based on the partial replacement 



design of Section 3. Assume (A1)-(A2)-(A3)-(A5). Then as N -4 oo, 



var (ise) ~ h IL liw^ 01 *'' 1 )' **■ 

Under additional assumptions, it is possible to find the asymptotic expression 
of Xh(t,t') in the previous theorem. Let G be an antiderivative of the density g 
in (A2). 

Corollary 1. Assume the conditions of Theorem^and (i) the sample sizes nh(-) 
are constant over time, (ii) limjv_ s>00 (a^m/fl — fh)) = < oo exists. Then as 
N ->• oo, 

E w l ~r^ exp ( ~ Ch " 7h(M,) ) ^ 
/i=i / 

The previous condition (ii) is reasonable since (a/ 1 m)/T is the average sample 
replacement rate per unit time, which in practice stays bounded. 

Under the assumptions of Theorems [3]|4] and Corollary [TJ we now compare 
the full and partial replacement designs in terms of variability of the ISE. As in 
Section 4 we include fixed panels as a special case of partial replacement where 
oth = Ch = 0. In comparison to fixed panels, partial replacement with Ch > 
induces an exponentially decreasing function in Var(ISE). The decrease is all the 
larger as the Ch are large and the data are (positively) correlated. In comparison 
to partial sample replacements, the full replacement design divides the order of 
Var(ISE) by a factor m, which massively stabilizes the estimation performance. 

Remark 1. If the survey objectives include the integral In = Jq HN(t)dt, then 
In = Jq P>ht(t)dt provides an unbiased estimator whose variance follows from 
the previous results. As before, in comparison to fixed panels, partial replacement 
of the sample reduces \&i(In) by an exponentially decreasing function and full 
replacement divides the order ofVar(lN) by a factor m. 

6. Composite estimation 

The estimation of ^iv(i) can be improved by using past data in addition to 
the current observations Xj-(t),k G s(t). Following the principle of composite 
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estimation (see e.g., Eckler, 1955), we utilize the partial replacement design of 
Section 3 and define a new estimator fi c (t) in an iterative fashion as a linear 
combination of the Horvitz-Thompson estimator jX^it), which is only based on 
the current observations, and of the estimator fi c (t — 5) updated with the esti- 
mated change in fj,jy between t — 5 and t, where 5 E (0,T) is a time window to 
be specified. Let < t < t' < T. If \v{t) — v{t')\ < 1, the estimator 



h=l VV H) n\ ) keSh{t)nSh{tl) / 

is unbiased for the population mean change A/iN(t,t') = ^v(i') — Atjv(*)- If 
\u(t) — v(t')\ > 1, the previous estimator is extended as 

Afj, N (t,t') = A[j, N (t,T u ( t ) +1 ) + An N (T r ,T r+1 ) + A^ A r(r i/ ( t , ) ,t / )- (6.2) 

r=i/(t)+l 

The composite estimator is then defined by fi c (t) = jlhtit) for t £ [0, ri) and 
Ac(t) = Q Aftt(i) + (i - Q) (Ac ((* - S)+) + A^((t - 5)+, t)) (6.3) 

for i G [t\, T], where Q £ [0, 1] must be specified and x + = max(x, 0). Note that 
if ah = with constant sample size rih(-) for all h, or if 5 = 0, or if Q = 1, then 
p- c (t) reduces to fcht(t)- I n this sense, the composite estimator can be thought 
of as a shrinkage estimator reverting to the Horvitz-Thompson estimator. The 
parameters ah, 5, and Q govern the tradeoff between past and present data in the 



estimation. The theoretical study of the estimator (6.3) and the development of 
automated selection procedures for the ah, 8 and Q are beyond the scope of this 
paper. However, important insights for selecting these parameters are provided 
in the next section. 

7. Numerical study 



This section examines the numerical performances of the HT estimator (2.2) 



and composite estimator (6.3) based on the full and partial replacement designs 
of Section 3. The simulations are based on real electricity consumption curves 
studied by Cardot and Josserand (2011). In this study the electricity usage of 
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N = 18902 French firms was recorded every half hour during two weeks. Based 
on the functional principal components analysis (FPCA) of the second week of 
data, we generate realistic new data as follows: 



3 

Xk(t) = f4f(t) + Zek + £k (*)> 
=i 



(7.1) 



where k G {1, . . . , N}, N = 10000 and t G [0, 1]. The term /j,* N in (7.1 ) is the mean 
function of the real data; the are independent random variables distributed 
as iV(0, <r|); o\,a\,o\ and 0i,</>2>03 are the first eigenvalues and eigenfunctions 



of the FPCA (see Figure 7.1). The process is a Gaussian white noise; it is 
independent of the Z^ and verifies Vax(sk(t)) = <5 2 Var(]T^ =1 Z^. <f>t(t)) with 
5 = 3%. The target function is [in = (1/N) X^fcLi -^k- 





Figure 7.1: Typical curves with mean function \x in thick solid line (left panel) and 
eigenfunctions <f>i,<f>2, ^3 (right panel) . 



In the simulations the curves are discretized at d = 400 equidistant 
points tj = j/(d + 1) and grouped in H = 5 strata according to their average 
levels Jq 1 Xk(t)dt. The strata cutoffs are determined by the box plot. We set the 
sample size to n(t) = 0.05 N at each instant £ £ [0,1] and consider two sample 
allocations: [PROP] allocation proportional to stratum size with nh(t) = nN^/N 
and [OPT] optimal allocation with n/ l (t) oc Nhjh(t, t) 1 / 2 (Neyman allocation). 



The strata sizes and sample sizes in each stratum are displayed in Table 7.1 (the 
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intervals in the bottom row are the range of nh(-) under optimal allocation). 



h 


1 


2 


3 


4 


5 


N h 


34 


2466 


5000 


2456 


44 


PROP: n h 


2 


123 


250 


123 


2 


OPT: n h (t) 


[1,2] 


[126,134] 


[226,240] 


[129,138] 


[1,2] 



Table 7.1: Strata sizes and sample sizes under different allocations. 

The replacement times are defined as r r = t^ r , r = l,...,m, with m = 
[d/3] = 133, which corresponds to a period of lhl5 between successive replace- 
ments. The sample s(-) is built either by full or partial replacement, with the re- 
placement rates a\ = ... = an '■= a and a G {0, 0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 1} 
under partial replacement. Each type of replacement (full or partial) is combined 
with each type of allocation (proportional or optimal). For each combination we 



approximate the mean and variance of the ISE for the HT estimator (2.2) based 
on 5,000 to 100,000 Monte Carlo simulations of s(-). 



Proportional allocation 



Optimal allocation 



Figure 7.2: Horvitz-Thompson estimator of the mean function. Mean value (solid line) 
and standard deviation (dashed line) of the ISE in function of the replacement rate a. 



In Table 7.1 the sample sizes under proportional and optimal allocation are 
quite close. This closeness stems from the data-generation mechanism (7.1) and 
the stratification. Indeed, in our simulations as in reality the main source of 
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variation for electricity usage is scale effects; in contrast longitudinal variations 
are qualitatively similar in the population. Also, the stratification used produces 
homogeneous strata with comparable standard deviations 7/i(i,t) 1//2 at a given 



time t. It is thus not surprising that in Figure 7.2 the estimation error is almost 
identical under proportional and optimal allocation. In the same figure, the MISE 
remains constant for all values of the replacement rate a as predicted by theory 
(see Section 4.3). A tiny reduction of the MISE can be achieved by using optimal 
allocation rather than proportional allocation; a much larger decrease would 
occur if the strata variances fluctuated differently over time. Importantly, Figure 
|7.2| illustrates the rapid decrease of Var(ISE) as a increases to 1 (see Corollary 
[TJ. In comparison to a fixed panel (a = 0), both the partial replacement design 
with a = 1 and the full replacement design reduce the standard deviation of the 
ISE by a factor 8. 

In the second part of the simulations we turn our attention to the composite 
estimator (6.3). Keeping the previous data, population stratification, sample 
sizes nh(t) and replacement times r r , r = l,...,m, we select the sample s(-) 
by partial replacement and proportional allocation. The replacement rate a 
varies in {i/10 : i = 1, . . . , 5} and the tuning parameters of the estimator satisfy 
5 G {1.25, 2.5, 6, 12, 24, 48} and Q G {i/10 : i = 0, . . . , 10}. Note that because the 
strata samples si(-) and s^(-) have size 2, a = 0.5 is the largest value for which 
(6.3) can be computed. To facilitate the interpretation, 5 is expressed in hours 
(recall that the survey duration is one week). 

The mean estimation error MISE is displayed in Figure 7.3 in terms of a, 5 
and Q. To simplify the visualization, two parameters vary in each graph with 
the third one held fixed. The figure reveals a strong increase in the MISE as 
Q — > 0. This increase is very fast for small S and slower as 5 increases. For 
all 5, values of Q close to 1 (that is, when the composite estimator is close to 
the HT estimator) yield reasonable performances. On the other hand, Table 



7.2 indicates that the optimal MISE is an increasing function of S (<5 = 25mn 



is the time step used for the discretization grid tj, j = l,...,d). This result is 
satisfactory in practice: it means that excellent performances can be obtained 
while retaining only very recent past data. This greatly alleviates data storage. 



Table 7.2 also shows that in order to obtain the best MISE, a must increase 
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and Q must decrease as 5 decreases. A possible explanation for the two previous 
findings is that as the composite estimator /t c (i) gets updated with respect to 
an increasingly distant time point t — S, the estimation of the population change 
A/jN(t — 5, t) becomes increasingly unreliable; to temper this, A//jv(* — S, t) must 
be based on more matched units (i.e. smaller a) and more emphasis must be 
placed on the estimation based on the current data (i.e. larger Q). The above 
analysis results for the MISE also hold for the variance of the ISE. 




Figure 7.3: MISE of the composite estimator in terms of a and Q with 5 fixed. 



s 


25mn 


50mn 


1.25h 


2.5h 


6h 


12h 


24h 


48h 


min QiQ (MISE) 


3.11 


3.26 


3.28 


3.60 


4.21 


5.10 


5.86 


6.21 




0.5 


0.5 


0.5 


0.5 


0.45 


0.3 


0.2 


0.05 


Qmin 


0.1 


0.2 


0.2 


0.4 


0.5 


0.6 


0.8 


0.7 



Table 7.2: Optimal MISE and arguments a, Q of the composite estimator in terms of 
S. 

In simulations not included here, the results of this section have been seen 
to hold with other types of stratification (e.g., fc-means) and other replacement 
frequencies m. Also, the full and partial replacement designs produce excellent 
performances for the estimation of the integral fj,pf(t)dt (see Remark lL 
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8. Discussion 

With this paper we have taken initial steps in the exploration of time- varying 
samples with functional data. We have developed rotation sampling designs as 
well as a composite estimation method. Their good statistical properties have 
been established in theory and in practice. In particular these methods have 
been seen to outperform fixed panels and the Horvitz-Thompson estimator. 

Hereafter we mention possible extensions of this work and related research 
directions. A first extension pertains to informative sampling. In addition to 
the adaptive rules mentioned in the introduction, an interesting possibility for 
the partial replacement design would be to let the replacement rates ah be ran- 
dom and time- varying. It would then be sensible to adjust oth(t) to the balance 
between transversal and longitudinal variations in the stratum Uh at time t: 
Oih(t) would be larger when transversal variations dominate longitudinal ones 
and smaller in the opposite case. A related problem is to find adaptive strategies 
that adjust the frequency of replacements at a given time (with our notations, 
the t t ) in function of the magnitude of transversal variations (see e.g., Marbini, 
2003). Of course, when developing a survey sampling strategy for a particu- 
lar application, the search for statistical optimality should be balanced with the 
various costs and constraints incurred by sampling. Regarding the composite 
estimator introduced in the paper, it would be useful to study its theory (in 
particular derive its variance) and to devise automated procedures for selecting 
its tuning parameters 5 and Q. Beyond composite estimation, effective survey 
estimation methods in the functional data setup could be based on functional 
linear models or nonparametric function regression approaches (see e.g., Ram- 
say and Silverman, 2005, and Ferraty and Vieu, 2006). Finally, our approach 
could be enhanced by incorporating auxiliary information in the survey. Model- 
assisted paradigms are available when the sampled curves are observed over the 
entire study period (Cardot, Chaouch, Goga, and Labruere, 2010) but methods 
for time-varying samples remain to be developed. In this regard it would be en- 
lightening to compare the benefits of using auxiliary information in design-based 
approaches (for survey weights) versus model-assisted approaches. 
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Appendix 

Proof of Proposition [2] 

The sum under study can be decomposed as Ylt=i Ae(t,t'), where 

A e (t, t')= E WMtWW)) XiWXiiftx^xtf) 

i,j,k,l£U h 

and Cijki = #{i,j,k,l}. To compute the Ag, we derive E (Ii(t)Ij(t)Ik(t')Ii(t')) 
based on the properties of SRSWOR and develop sums ^ Xi{t)Xj(t)X^{t')Xi{t') 
using the identity Ylk^u h Xk(t) = 0- Let i*,j*,k*,l* be four distinct units in Uh- 
We begin with the straightforward calculation of Ai(t,t'): 



(1) 



The term A%(t, t') can be expressed as follows: 



A 2 {t,t') =E (/*.(*) Jfc.(0) 



(N h -l) 2 lh (t,t) lh (t',t' 



Y X 2 k (t)Xl{t') 

k&U h 



2E(J i .(t)Ii*(t / )4*(t)4*(t')) 



k&U h 



- 2 [E (I*. (4)4.(0^.(0) + E(l i *(t)I k *(t)I k *(t')) ] E **(*)**(0- 

(2) 

Next, it can be shown that 

A 3 (M') = [E(l f . (*)/,-. (tJiHO) + E(l i .(t)/ fc .(0ij.(0)] 

x -( Nh -l) 2 lh (t,th h (t',t') + 2 Y X 2 k (t)X 2 k (t')] 
+ 4E(Ji.(*)J i .(0/i.(t)/*.(0)[-(JV fc -l) 2 7g(t,0 + 2 E **(*)**(0 
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To compute A±(t,i?), use the decomposition 



J2 x i (t)x J (t)x k (t')x l (t') = Yl Yl M^Xj^Xk^Mt') 



i,j,k,leU h 



=1 i,j,k,leU h 



together with ([TJhQ to obtain 

A 4 (t,t') = E(l i * (t)Ir (i)4*(i'K**(*')) x 



{N h -l) 2 lh {t,t) lh {t',t') + 2{N h -l) 2 1 2 h {t,t')-Q Y, Xk(t)Xk(t') 

keu h 



(4) 



The proof is completed by gathering Q-Q and observing that all terms 
involving ^keu h ^fcW^f (0 are of negligible order 0(N h ) thanks to (Al). □ 

Proof of Proposition [3] 

Let k,l be two distinct units in a stratum Uh- Consider the Markov chain 
{(Jjt + J;)(r r ), r = 0, . . . , m} which counts how many units among k, I are present 
in the sample at the successive replacement times. This chain has three possible 
states: 0, 1, and 2. For r = 1, . . . , m, the transition probability matrix 



(P + J,)(T P ) = j - 1|(4 + WiTr-l) 



1)) 



1<*J<3 



can be represented as 



where 



T3* 



p; + E r , 



2 (1 - fa) f3 r 



(5) 



a & (l-A-) aifcA- + (1 - a h ) (1 - A-) (l-a h )A 



V 



or 



2 (1 - a/j) 



(1 - a?t) z 



/ 



and p r = P (k £ s h {T r )\k £ s h (r r -i)) = (A(r P ) - (1 - a h )/fc(Tr_i)) / (1 - /ft(r r _i)). 
Recall that = P (A; ^ s/i(T r )|A; 6 s^(r r _i)). The matrix E r , whose cumbersome 
expression is not given here, is asymptotically negligible in comparison to P*. 
More precisely, (A5) guarantees that max r= i v .. im ||E r || = O (1/iV/j) as JV -> oo, 
where || • || denotes an arbitrary matrix norm. For simplicity we use the spectral 

1 /2 

norm ||A|| = sup x ^ 

(x'A'Ax/x'x) 1 '' henceforth. 
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The transition probability matrices P r have unit spectral norm. Using the 
binomial formula, the triangle inequality, and the inequality ||AB|| < ||A|| • ||B|| 
holding for all compatible matrices A and B, it follows that 



v(t>) v{t>) 

n p .- n 

r=u(t)+l r=u(t)+l 



II Pr- n P * ^ E f^^-^)- 1 )/ max ||E 

< ( 1 + max ||E r 
r=l,...,m 



= 0\m max ||E r || (6) 

V r=l,...,m J 

uniformly in < t < t' < T. Combining the previous results and condition (A5), 
it comes that 

v{t') «/(*') 

n p r = (i+ (i)) n p * ( 7 ) 

r=v(t) + l r=v(t)+l 

We now study the simpler product E[r=i(t)+i P *- Without loss of generality, 
set v(t) = and v(t') = r. Define Q r = IT[=i an d write [Mij f° r tne j)th 
coefficient of a matrix A. It remains to compute [Q r ]i3, [Qr]23> and [Q r ]33- 

Using the facts that Q r is a transition probability matrix and that [P*] 2 j = 
[P*]^ 2 [P*]^ 2 for j = 1,3, it can be shown by induction that [Qr]ii 2 + [Qr]i3 2 = 1 

1/2 1/2 

for r = 1, . . . , m. Consequently [Q r ] 13 = (1 — j3 r — a^) [Q r -i]i3 + which can 
be reformulated as 

[Qr]i3 2 " AM = l : ah - fh{T : ] ([Qr-tii 2 ~ h{r r A . (8) 

1 - Jh\Jr-l) V / 

1/2 

Noting that [Qi] 13 - fain) = -fh(T )(l-a h -fh(n))/{l-fh{To)) and iterating 
we obtain the identity [Qr]}^ = fh(r r ) - fh{ro) A/ l (r , r r ). 

1/2 1/2 

The arguments used to determine [Q r ]i3 can be identically applied to [Qj+33 • 

1/2 

Omitting the calculations, we directly state that [Q r ]33 = (1 — /fe( r o)) -Wto, r r)+ 
fh{T r )- The expressions of [Q r ]i3 and [Qr]33 can be checked by induction. 
Finally we turn to [Qr]23- The total probability formula yields 

F(k,l G s(r r )) = F(k,l G s(r r )\k,l G s(0))F(k,l G s(0)) 

+ 2F(k,l G s(T r )\k G s(0), / i s(0))F{k G s(0), I £ s(0)) (9) 
+ P (k, I e h s h {t')\k, I i s(0)) P (fc, Z i s(0)) . 
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In view of ([7]) we obtain the approximation 

fh(r r ) 2 ~ (A(t ) 2 [Qr] u +2 f h (r ) (1 - f h (r )) [Q r ] 21 + (1 - A(r )) 2 [Q P ] 31 ) . 

(10) 

The proof is completed by plugging the expressions of [Q r ]i3 and [Q r ]33 in 
the previous relation. □ 

Proof of Theorem [3] 

We start by finding the asymptotic expressions of Ci(t,t') and C2(t,t') in 
Proposition [2| Exploiting (A3), the basic properties of SRSWOR and the inde- 
pendence of the s/j(r r ), r = 1, . . . , m, under the full replacement design, it comes 
that 

f C X (M') ~ hit) f h {t') (1 - Ut)) (1 - hit')) 

\ C 2 {tJ) ~2f h {t)f h {H) (1 - f h (t)) (1 - hit')) 8 v{t)v{tl) 

uniformly in t,t' £ [0, T] as N — > oo. Hence the second term in the right-hand 
side of (5.3) cancels out with the term in Ci(i, if) of Proposition [2j 
Writing Var(ISE) = (2/7V 2 ) JL T]2 cp N it,t')dtdt' , we deduce that 

uniformly in t, t? 6 [0, T] as JV oo. Hence, using (A2), (A4) and the mean 
value theorem, we get 

2 x - N h N h , ^ (1 - hjr r )) (1 ~ Mtv)) /•/• ,,,,,, 
Var(ISE) ~ ^ g g ^ ^ t ) 7 ^, * )^ 

2 ^ N h N h , ^ (1 - hjr r )) (1 ~ Mtv)) 2 



^ 2 ^ ^ 2 ^ A(r.)/^(r r ) 

2 ^ N h N h , ^ (1 - hjrr)) (1 ~ MtvQ) g(r r ) 2 
^ ^ ^ ^ hirr)h'irr) m 2 

S^io /UW) ^)7,(M)7.(M)^. □ 



mJV 2 ^ iV 2 
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Proof of Theorem |4] 

This result is established along the same lines as Theorem [3| We start by 
finding the asymptotic expression of C\{t,t') and C2(t,t') in Proposition [2] under 
partial replacement. In view of Lemmas [Tj(2] and Proposition [3j it holds that for 
all < t < t' < T, as N -> 00, 

E(7 i .(t)/ fc .(0) 

= P (k* G sCOl**,*!* G s(t))P(r,fc* G s(t))+P(fc* G s(i')|«* G s(t),fc* $ s(i))P(i* G s(i), F g «(t)) 
= P (fc* G s(t')\k* G a(t)) P («*, k* G s(i)) + P (k* G a(0|fc* £ s{t)) P (i* G s(f), k* <£ s(t)) 

~ [(i-A(t))A h (t,0 + A(0 ]/£(*) + [A(0- AWAh(t,t')]A(i)(i-AW) 
= h(t)h(t'). 

By symmetry this expression holds for all G [0, T]. Similarly, we find that 

E(/ i .(t)/ i .(t)/ fc .(0) ~ ft(t)f h (t), 

E(I t *(t)I k *(t')Ii*(t')) ~ fhWfttf), 

E(I t *(t)I r (t)I k *(t')Ii*(t')) ~ fl(t)f 2 h (t'), 

E (J*. (t)7fc, (f)) ~ [ (1 - f h (t))X h (t, t') + / fc (f ) ] 2 4 2 W, 

E(7 i »(t)7 i *(t')7 i «(t)7 fc .(t')) ~ [(l-A(t))A h (t,0 + A(0]4 2 (i)A(0- 
Therefore 

f C X (M') ~ / ft (t)/ h (f) (1 - (1 - 

\ C 2 (t,0~2/2(t)(l-/,(t)) 2 A2(t,t') 
Combining the previous result and Proposition^ it stems from (5.3) that 

]T ^(ttfXiiftXjiftXkWxtf) ~ (1 - f h {t)f\l{tM{t,t')Nl. 

i,j,k,leU h 

(12) 



On the other hand the inter-strata contribution to (5.2) can be simplified 
thanks to Theorem [2} 

, £ ») = N » W °' (13) 

With the notations of the proof of Theorem [3j one deduces from (12) and 
@ that 

Ij^^p^'^^'^J (14) 
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for all t,t' £ [0, T], as N — > oo. To apply the dominated convergence theorem, it 
suffices to check that the 4>n, N > 1, are uniformly bounded on [0, T] 2 . In view 



of (Al) and (A2), the right-handside of (14) has a finite number of terms, the 
terms (1 — fh(t))/fh(t') and 7^(4, t') are uniformly bounded with respect to h, 
t,t' G [0, T] and N, and \Xh{t,t')\ < 1 as a product of eigenvalues of transition 
probability matrices. The dominated convergence theorem thus applies, which 
concludes the proof of Theorem [4j □ 

Proof of Corollary [l] 

In view of Theorem|4j it suffices to show that \h{t, t') ~ exp (— \ G(t) — G(t')\) 
as N — > oo. By condition (i) of the corollary, we first see that A/ l (t, tf) is equal to 
(1 - a h /(l - fh)) luit) ~ v{t ' )l - By condition (A2), the term \v(t) - v(t')\ works out 
as (m + 1) \G(t u ^) — G(T v / t n\, which in turn is asymptotic to m\G(t) — G(t ; )\. 
Finally exploiting condition (ii) of the corollary and the approximation ln(l— x) ~ 
(— x) as x — > 0, we obtain the sought equivalent for Xh(t, tf). □ 
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